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ABSTRACT 

We present numerical investigations into the formation of massive stars from turbulent 



cores of density gradient p oc 



-1.5 



. The results of five hydrodynamical simulations 



are described, following the collapse of the core, fragmentation and the formation of 
small clusters of protostars. We generate two different initial turbulent velocity fields 
corresponding to power-law spectra P oc k~ 4 and P oc /c -3 - 5 , and apply two different 
initial core radii. Calculations are included for both completely isothermal collapse, 
and a non-isothermal equation of state above a critical density (10 _14 gcm -3 ). Our 
calculations reveal the preference of fragmentation over monolithic star formation in 
turbulent cores. Fragmentation was prevalent in all the isothermal cases. Although 
disc fragmentation was largely suppressed in the non-isothermal runs due to the small 
dynamic range between the initial density and the critical density, our results show 
that some fragmentation still persisted. This is inconsistent with previous suggestions 
that turbulent cores result in the formation of a single massive star. We conclude that 
turbulence cannot be measured as an isotropic pressure term. 
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1 INTRODUCTION 

There are several potential difficulties in forming high-mass 
stars. Firstly, the timescale of less than 10 6 years to 
assemb le 10 to more than 100 M implies large ac c retion 
rates (IZinnecker et ah| Il993l : iBehrend fc MaedeJ l200ll: 
iNorberg MaedeJ 1200(1) . Secondly, their c rowded loca- 
tion in the centre of young stel la r clus te rs (IClarke et al 
2000; Hillenbrand &: Hartmannl Il998t ICarpenter et al 
1997]; lLada Ladal 1200.^ limits the final mass of 
any collapsing fragment: confining the Jeans radius 
to be less than the interstellar separation means 
that the resultant Jeans mass will be correspond- 



ingly small (|Zinnecker 



Bonnell. Bate & Zinnecker 1998). Lastly, and most im- 
portantly, the radiation pressure from a high-mass star is 
sufficient to rev erse the infall of gas that contains typical 
dust properties (|Wolfire fc C assinelli 19871: iBeech &; Mitalaa 
1994). There are a number of ways this last problem can be 
circumvented. The radiation p ressure could be overw helmed 
by ultra-high accretion rates (McK ee &: TanlEo03l) . Alter- 
natively, accretion could occur preferentially through an 
equatorial disc, while the central object is rapidly rotating 
and thus emitting most o f its radiation towards the poles 
JYorke fc Sonn halter 2002h . A third solution is that massive 
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stars for m due to stellar mergers in the ultra-dense core of 
a cluster (Bonnell et alJll99fl Bonnell &: Batell2002h . 

The lMcKee fc Tanl fcOOlfeOOfl model depends on high 
accretion rates that would be expected to occur in a dense 
gas core in the centre of a cluster. This core is envisioned 
to be supported by turbulence, as thermal pressure is inade- 
quate, and adopts a steep density profile in order to prevent 
fragmentation. Even neglecting how such a core could arise 
in the centre of a stellar cluster, one potential difficulty is 
that turbulent clouds are known for their tendency to frag- 
ment and form a stellar cluster. It is this possibility that we 
address here. 

Fragmentation during the initial stages of star 
formation has been illustrated through numerical 
studies of collapsing turb ulent molecular clouds. 
Bate, Bonnell &; Bromml (|2003f) describe the collapse 
of a uniform 50M© turbulent cloud of diameter 0.375pc, 
which forms 3 cores containing 23+ protostars and 
27+ brown dwarfs. Fragmentation has also been 
demonstrated in isolated, rotatin g 1M cloud cores 
(Burkert, Bat e fc Bodenheimerl Fl997V where low multiple 
systems are produced despite an initial p oc r _1 density 
distribution. 

T he formation of stellar clust ers from turbulent cloud 
cores (IBonnell. Bate fc Vind l2003 ) has shown how the frag- 
mentation of the turbulently induced filamentary structure 
produces hundreds of stars. These stars fall into local poten- 
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tial wells forming small sub-clusters which eventually merge 
to form a large stellar cluster. The massive stars form in the 
centre of the sub-clusters due to competitive accretion, nat- 
urally explaining why mass ive stars are found in the cores 
of dense stellar clusters (iBonnell. Vine fc Batd 12004 V The 
stellar densities in these cores can b e very high, approach - 
ing that for stellar mergers to occur feonnell fc Rate 2002). 

Our calculations adopt a turbulent core with a density 
profile and velocity dispersion comparable to lMcKee &; Tar] 
( 2002j) . The purpose of this paper is to investigate whether 
the steep density gradient assumed by McKee and Tan is 
sufficient to prevent fragmentation and consequently provide 
a suitable approach to massive star formation. 



2 COMPUTATIONAL METHOD 
2.1 SPH code 

We use the 3D smoothed partic le hydrody namics (SPH) 
code based on the version by Benz (|Benzll990i) . The smooth- 
ing lengths between particles are allowed to vary, but the 
typical number of neighbours for each particle is N ne i g h ~ 
50. Artificial viscosity is included with the standard parame- 
ters a = 1 and [3 — 2. Gravitational forces are calculated us- 
ing a binary tree. The code uses 10 6 particles (so for a 3OM 
core (Section 2.3) the minimum particle mass is 3x 10 _5 Mq) 
and simulations were run for approximately half a free-fall 
time. All computations were performed using the United 
Kingdom's Astrophysical Fluids Facility (UKAFF), a 128 
CPU SGI Origin 3000 supercomputer. 



2.2 Sink particles 

Simulations in SPH become very computationally expensive 
as the density increases du ring the collapse of a core. T he 
insertion of sink particles (|Bate. Bonnell fc P rice 1995|) at 
a certain density is widely used to extend calculations. To 
resolve the local Jeans mass requires ~ 2N ne i g h particles 
feate fc BurkertJll99/t) . i.e. 3 x 10~ 3 M o in these simula- 
tions. This is then the minimum resolvable mass that can 
be considered unstable for collapse. The minimum density 
for insertion of a sink particle can then be determined from 
Jean's equation, inserting 3 x 10 _3 Mq as the Jean's mass: 



Mj 



5R g T \ 
2pG J 



3/2 



-1/2 



= 3 x 10~ d M ( 







(i) 



Taking the temperature T to be 20K and p — 2.46, this rear- 
ranges to give pace — 1.06 x 10 -13 gcm -3 . The corresponding 
radius for accretion IS T 'acc ^ 20AU. 



2.3 Initial conditions 

McKee & Tan (hereafter MT(2003)) assume a centrally 
condensed core, following a density profile of p oc r~ kp . 
We take k p = 1.5, the fiducial value of MT(2003). This 
is shallowe r than th e singular isothermal sphere model 
(k p = 2, (|Shul llQfjfO but comparable with observations 
of cloud cores (iMvers Fulled (|l992h give k p = 1.6, 



oi cloud cores (IMyer sj^buJ 
IWard-Thompson et al.| 



kn 



1.2). 



0.06pc with a mass of 3OM . With a resolution of 10 6 parti- 
cles, the initial density profile p oc r _3//2 rendered a central 
density of ~ 3 x 10 -15 gcm -3 , whilst the average density was 
pav ~ 2.5 x 10 -18 gcm -3 . The results are given in terms of 
the free fall time, tff = (3tt /32Gp av ) ~ 4.5 x 10 4 years. 
MT(2003) use similar initial conditions, but take a mass of 
60M© to form a 3OM star with 50% efficiency. We also 
performed computations wi th a core radius .2pc, since this 
gave similar dimensions to iBate et alJ (|2003). The average 
density was then 3 x 10 -19 gcm -3 and the central density 
2 x 10 -16 gcm -3 . The free fall time for these results was 
2.7 x 10 5 years. In all calculations, the initial temperature of 
the core was 20K, comparable to the cold gas of which the 
first massive stars form. 

As in equation (1), the Jeans mass is 



Mj(p(r)) 



5R g T 



3/2 



-7Tp 



-1/2 



(2) 



2pG J V 3 

where p the mean molecular weight is taken as 2.46. The 
number of Jeans masses contained in a core is then 

rR p(r)Anr 2 



Nj{R) = I 
Jo 



Mj(p(r))— (3) 

This gives 69 Jean's masses when R — 0.06pc and 11.5 when 
R 0.2pc. 

2.3.1 Turbulence 

MT(2003) incorporate turbulence by means of an effective 
turbulent pressure term whereby P oc r~ kp . By assuming 
hydrostatic equilibrium, kp = 2(k p — 1), and the resulting 
velocity dispersion relation is derived as 



a (X. c s — 



r (2-fc P )/2 



(4) 



where a is the velocity dispersion, and c s the effective 
isothermal sound speed. Inserting k p — 1.5 then gives 
a oc r 0-25 . This is shallower than observational results ob- 
tained for the similar Larson relation a oc L a (La typical 
length sc ale), which in dicate th at 0-25 < a < 0.75, e.g. 
a 0.38 (lLarsonlll98Jl . a 0.5 (|Mverdll983h . 

We simulate turbule nce using a grid method outlined 
by iDubinski et alJ ()l995l) and briefly described here. The 
initial velocity field is calculated from a random Gaussian 
distribution and follows a power spectrum 



P{k) = < \v k \ > 2 oc k~ 



(5) 



where k is the wavenumber. To obtain a divergence free ve- 
locity field, velocities are constructed from a vector poten- 
tial A, such that v = V x A. We sample components of 
Ak, the Fourier transform of A, from a grid of co-ordinates 
(fee, k y ,k z ). In order to generate a Gaussian distribution sat- 
isfying (5), components of Ak are selected as 

A k (fc, k y ,k z ) = |k| ^ (C kx e i9k * , C ky e^y , C ha e i0 ** ) 

(6) 

The amplitudes C are sampled from a Rayleigh distribution 
and phases < < 2tt from a random distribution (see 
IDubinski et alJ (^995)). Taking the inverse Fourier trans- 
form of V x Ak then gives 



The first computations took the core radius R to be 



2tt 



x A 



k e 



(7) 
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Numerically the real x component of (7) becomes 



3 RESULTS 



1 kmaX -n-2 

2^ I] \ k \~^ [k z C ky 8m(k.r + (f> ky ) 



-kmax 

■ k y C kz sin(k.r + (j) kz 



(8) 



summed over Kx i k y and k z . Similar expressions provide v y 
and v z , for a a particle at position (x,y,z). The parameter 
kmax in equation (8) is often chosen to satisfy the resolu- 
tion requirements of the model, i.e kmax ~ (2tt / L m in) where 
Lmin is the smallest length scale or minimum separation be- 
tween particles. However, since our model is centrally con- 
densed and the density varies considerably, we allow kmax 
to vary for each particle according to the local smoothing 
length. We take 



{kn 



(9) 



hi is the smoothing length of the zth particle. In this way, 
velocities in less dense (low resolution) regions, could be 
computed much faster than those in the dense centre of the 
core (high resolution), reducing the computational expense 
of the calculation. 

The velocity field generated gives a dispersion - length 
scale dependence similar to the Larson relation. Calcula- 
tions are included for a oc r 0-25 , consistent with MT(2003), 
and a oc r 0-5 which corresponds better with observations. 
These dispersion laws are attained by generating power 
spectra with n = 3.5 and n — 4 respectively. However 
iMvers &: Gamrniel (|l999T) show that n ~ 3.25 provides a bet- 
ter fit for the a oc r 0-25 dispersion law. 

The core is initially in virial equilibrium, such that 

Egrav + 2Eturb + ^Etherm = 0. (10) 

The ratios of turbulent to gravitational energy were 

Eturb 



J-^grav 
Eturb 



0.4 for R = 0.06pc, 



0.25 for R = 0.2pc. 



(11) 
(12) 



The turbulent energy is then all owed to dissipate ove r the 
dynamical timescale of the core (iMac Low et al.lll99/fl . We 
do not include any driving, which if included, should increase 
the susceptibility to fragmentation. 

2.3.2 Equation of State 

The thermal behaviour of the cloud core was assigned as 
either: 

(i) Completely isothermal. 

(ii) Non- isothermal above a critical density, allowing the 
gas to heat up as the core collapses. In this case the equation 
of state exponent varies as 



7 



1 if p < 10" 14 gcm- 
1.67 if p^ 10 _14 gcnT 



(13) 



(iii) Non- isothermal above the critical density with a sec- 
ond isothermal collapse phase for p > 10 -12 gcm -3 . 

Although the core is initially globally in virial equilibrium, 
the inner regions are actually over-supported due to the 
isothermal nature of the gas. 



A summary of the calculations undertaken, their parameters 
and the overall results are shown in Table 1. The parameter 
a describes the initial velocity dispersion (where a oc r a ). 

In all simulations, turbulence initially supports the core 
in virial equilibrium. As the simulation progresses, turbu- 
lence decays through shocks in the core, dissipating kinetic 
energy. The core then begins to undergo collapse, with den- 
sity increasing until the formation of individual protostars. 
The morphology of the cloud is similar in each simulation - 
the core becomes elongated, generating significant structure 
which forms the basis for subsequent fragmentation (e.g. 
Figure 1). 



3.1 Isothermal results 

Models 1, 2 and 3 used an isothermal equation of state 
throughout the simulation. Models 1 & 3 apply a a oc L 0-25 
dispersion law whereas Model 2 uses a oc L ' 5 although the 
initial velocity fields are determined by the same series of 
random numbers in all calculations. Model 3 differs by tak- 
ing a 0.2pc radius cloud. 

All 3 simulations produced significant numbers of pro- 
tostars. Table 1 includes the total number of protostars 
formed, although results cannot be directly compared as 
some formation was still ongoing. 



3.1.1 Comparison between dispersion laws 

Figure 1 displays different stages during Models 1 & 2, 
which apply different power spectra. The main difference be- 
tween Models 1 & 2 is the timescale for protostar formation 
and fragmentation. Collapse and protostar formation occur 
later in Model 1 compared to Model 2, although their over- 
all evolution is similar. This is due to the shallower power 
law P oc A: -3-5 (compared to P oc k~ A for Model 2 when 
a — 0.5) which supplies more kinetic energy over smaller 
length scales. This provides greater support over small scales 
and generates more structure. Consequently density plots for 
Model 1 show more structure than those for Model 2, and 
the protostars are less widely distributed (Figure 1). 

The left side of Figure 1 illustrates the evolution for the 
a = 0.25 dispersion law, which best represents MT(2003). 
A protostar forms in the central part of the core after 
0.19£//, surrounded by a disc of gas. The core fragments 
after 0.25£//, when 2 further protostars form from the disc 
material. A fourth protostar is also created via fragmenta- 
tion of an elongated filament extending from the initial pro- 
tostar (Figure 1 ii). No further protostars form for ~ 0.15£//, 
as the core is largely supported by internal rotation. How- 
ever the internal dynamics are still chaotic and further col- 
lapse occurs as support is lost. Further fragmentation occurs 
through collapse along filaments near the core's centre (Fi- 
grure 1 iii). This results in a second phase of star formation 
at around 0.47£//, producing a total of 28 protostars after 
0.54*/,. 

Model 2 (Figure 1, right), where a = 0.5, shows similar 
behaviour. Again an initial protostar is formed after 0.09£// 
accompanied by a disc which subsequently distorts and frag- 
ments, producing 7 protostars after 0.13£//. Similarly there 
is a second phase of fragmentation through the collapse of 
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i)t=0,24 ii : 0.11 




Figure 1. Column density plots (logarithmic scale, minimum and maximum densities of 4 x 10 _17 gcm 3 and 4 x 10 _ 13 gem 3 ) showing 
evolution of central region of the core for a oc L 0,25 (left) and a oc L 0,5 (right). Length scale for each plot is 0.02pc X 0.02pc. The number 
of protostars in the plots are: Left - i)l, ii)4, iii)6, iv)13; Right - i)l, ii) T, iii)7, iv)14 . Times are given in units of tff = 4.5 x 10 4 years. 
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Model 


a 


R (pc) 


Jean 


Predominant time of 


No. of 


Total mass 


Most massive 








No. 


formation (£//) 


protostars 


accreted (Mq) 


protostar (Mq) 


1 


0.25 


0.06 


95 


0.47 


28 


2.73 


0.435 


2 


0.5 


0.06 


95 


0.37 


14 


2.12 


0.528 


3 


0.25 


0.2 


16 


0.55 


19 


1.79 


0.66 


4 


0.25 


0.06 


95 


0.54 


3t 


2.02 


0.99 


5 


0.25 


0.2 


16 


0.54 


2 


2.10 


1.28 



Table 1. Table showing results of all simulations, giving the total number of stars formed and the approximate time at which most 
formation occurred. Models 4 Sz 5 use a non-isothermal equation of state above a critical density (Section 2.2.2). However Model 4 also 
includes a second isothermal collapse phase. ^ Further fragmentation apparent but sink formation suppressed by equation of state. 



elongated filaments, leading to a total of 14 protostars after 
0.44£//. Figure 2 plots the mass accreted and the number of 
protostars formed against fraction of a free fall time, for both 
dispersion laws. Models 1 and 2 display very similar profiles, 
offset by ~ O.lt//. The number of protostars formed in each 
case shows 2 periods of star formation at 0.1/0.2t// and 
0.38/0.48t// separated by a period of 0.26/0.15t//. 

3.1.2 Overall comparison of isothermal runs 

The isothermal results all display a high level of fragmenta- 
tion, demonstrating that turbulent pressure is anisotropic. 
Turbulence cannot be represented by an isotropic thermal 
pressure, instead distortion of the core leads to the formation 
of multiple protostars. The degree and timescale of fragmen- 
tation vary according to the turbulent power law and the 
size of the core. Whereas a steeper velocity-sizescale rela- 
tion (a — 0.25) advances protostellar formation, collapse is 
delayed for the R = 0.2pc core in Model 3. The larger core 
has more thermal support so the ratio of thermal to grav- 
itational energy is greater (equations (11), (12)). Protostar 
formation begins at ~ 0.52£//, so collapse of the core and 
subsequent fragmentation occurs at a higher fraction of the 
free fall time compared with Model 1 where R = 0.06pc. 
Correspondingly, density plots for Model 3 also show less 
structure, tending to retain a centrally condensed profile, at 
least over large length scales. There was, however, still suffi- 
cient structure over the centre to produce 19 protostars over 
0.6t//, so significant fragmentation still occurred. 

3.1.3 Accretion rates 

The main motivation for MT(2003) was to achieve a high ac- 
cretion rate to overcome radiation pressure. Accretion rates 
can generally be estimated as the resulting protostellar mass 
divided by the free-fall time (M/t ff ). Thus MT(2003) con- 
struct their models to produce an early accretion rate of 
lO~ 4 M 2/r _1 , which increases to 10~ 3 Moyr' 1 . From Fig- 
ure 2 we see the accretion rate is initially ~ 10 -4 M^yr -1 in- 
creasing to 10 -3 M02/r _1 after 0.3£//. These accretion rates 
agree with MT(2003), but are based on the total mass ac- 
creted onto many protostars (between ~5-20). The accretion 
rate for an individual protostar is therefore significantly less. 
In contr ast, the cluster accre tion models for massive star for- 
mation (iBonnell et all 12004 ) show that although the mean 
accretion rate (M/tff) is only 10 _6 M©?/r _1 , the actual ac- 
cretion rate onto the growing massive star in the centre of 
a cluster is 10 -4 M©?/r _1 . 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 
Timeft ) 



Figure 2. Plots for number of stars (top) and mass accreted (bot- 
tom) versus time for Model 1, a = 0.25 (solid line) and Model 2, 
a = 0.5 (dashed line). 



3.2 Non isothermal equation of state 

The large scale collapse of the non-isothermal models are 
very similar to the corresponding isothermal cases, since 
the equation of state changes only in the central denser 
regions. However, Table 1 shows that the number of pro- 
tostars formed is dramatically reduced (Models 4 & 5). In 
Model 4, (0.06pc core), the equation of state changed be- 
fore fragments could collapse to form sink particles. A sec- 
ond isothermal collapse phase subsequently produced 3 pro- 
tostars. For Model 5 (0.2pc core), the core is initially less 
dense so 2 protostars were able to form without any second 
collapse phase. 
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To understand the difference between protostar forma- 
tion in the isothermal and non-isothermal models, it is nec- 
essary to compare the protostar-forming region in each case. 
Figure 3 illustrates the contrast between Model 1 (isother- 
mal) and Model 4 (non-isothermal) in the dense part of the 
core. Complex structures continue to exist on small scales in 
the isothermal case, including a disc and close companions 
formed through disc fragmentation. In comparison, Model 4 
is dominated by 2 large volumes of dense thermally sup- 
ported gas. 

On scales comparable to the core radius however, den- 
sity profiles show that considerable structure still exists in 
the non- isothermal core. In Figure 4 (Model 4), the maxi- 
mum densities over the x and y directions are plotted (while 
the equation of state remains non- isothermal). Five peaks 
have been selected from the profiles which show a large den- 
sity contrast with the surrounding material. These indicate 
independent, disconnected fragments. Figure 5 shows an ac- 
companying column density plot, obtained at the same time 
frame, which illustrates the large scale structure of the core. 
The locations of the five peaks are indicated, showing where 
fragmentation is taking place. The 2 highest peaks in the 
densest part of the core correspond to the dense regions 
shown in Figure 3. These also subsequently form 3 proto- 
stars during the second isothermal collapse phase, as the 
larger region (of Figure 3) sub- fragments. With a longer pe- 
riod of time, and a less severe equation of state, we would 
expect all the dense fragments corresponding to these peaks 
to collapse independently and form protostars. 

The distortion of the 0.06pc core suggests that mono- 
lithic formation is unlikely, despite the formation of very 
few protostars. Our analysis suggests fragmentation will oc- 
cur producing at least 6 separate bodies (5 peaks identified 
from Figures 4 & 5, and isothermal sub- fragmentation of 
peak 2). This is again a consequence of anisotropic turbu- 
lence dominating homogeneous thermal support. 

Conversely, the larger core (Model 5) remains predom- 
inantly centrally condensed. Thermal pressure supports the 
core over large scales, whilst the equation of state prevents 
any local fragmentation. With only 2 centrally located pro- 
tostars, Model 5 would only be expected to produce a low 
multiple (2 or 3 protostars) system. 



4 CONCLUSIONS 

Turbulently supported clouds of typical density structure 
p oc r -1 ' 5 have been hypothesised to be the progenitors 
of massive stars. Their density concentration could pre- 
vent fragmentation and their high accretion rates over- 
whel m the radiation pressure from the accreting massive 
star JMcKee fc Ta n 2003, 200^). We have performed numer- 
ical simulations of this scenario and found that the initial 
centrally condensed density profile proved insufficient to pre- 
vent fragmentation. Turbulent support generates significant 
structure in the core which forms the basis for subsequent 
fragmentation. Turbulence cannot be assumed to act as the 
equivalent of an isotropic pressure. In addition, although the 
total mass accretion rates are comparable with MT(2003), 
individual protostellar accretion rates are significantly lower. 

The fragmentation of the centrally condensed 30Mq 
cores forms ~ 20 stars over the time period of our isother- 




Figure 3. Logarithmic column density plots (as Figure 1, min- 
imum and maximum densities of 4 x 10 — 17 gcm — 3 and 4 x 
10 -i3 gcm -3) when a = o 25, R=0.06pc for isothermal (Model 1, 
top) and non-isothermal (Model 4, bottom) cases. The size of the 
plots are 0.02pc by 0.02pc and the time is 0.32t //; t ff = 4.5 X 10 4 
years. 



mal simulations. The protostars form from a combination of 
the fragmentation of filamentary structure and the fragmen- 
tation of rotationally supported discs. Our non-isothermal 
runs suppress this disc fragmentation as the gas is assumed 
to heat on scales of the discs. Nevertheless, additional frag- 
mentation is expected once a second collapse phase softens 
the equation of state. Competitive accretion i n these clusters 
will determine the eventual stellar masses (Bon nell et alJ 
I2001T) . 

These simulations neglect the potential effects of any 
magnetic fields present in the core. We do note that this 
should not impede the fragmentation process as magnetic 
fields have been shown not to affect t he structure gen- 
eration in turbule nt molecular clouds ([Stone et al.l 1998; 
IPadoan et all2004h . Even in the presence of magnetic fields, 
turbulence does not act as an isotropic support, and neither 
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-4 -2 2 A 

x,y, (0.01 pc) 

Figure 4. Logarithmic plot of maximum density (over all points in the simulation) versus x (red) and versus y (green) for Model 4. 
Dotted lines show 5 dominant peaks, p measured in computational units where 1 unit=7 x 10 _17 gcm -3 . The plot is taken at t = 0.53tff 
(tff = 4.5 x 10 4 years) when the equation of state is non-isothermal. 




Figure 5. Column density plot for Model 4, with peaks identified 
in Figure 4 indicated. Only points where log(p) > —2 have been 
selected to emphasise the structure. The size scale of this plot 
is O.lpc by O.lpc with minimum and maximum densities of 7 x 
10 -19 gcm -3 and 4 x 10 -13 gcm -3 . Again the plot corresponds to 
a time of t = 0.53tff. 

turbulence nor magnetic fields should be modelled as an 
isotropic pressure term. 
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